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SUMMARY 

A multigrid preconditioned conjugate gradient method (MGCG method), which uses the 
multigrid method as a preconditioner of the PCG method, is proposed. The multigrid method has 
inherent high parallelism and improves convergence of long wavelength components, which is 
important in iterative methods. By using this method as a preconditioner of the PCG method, an 
efficient method with high parallelism and fast convergence is obtained. First, it is considered a 
necessary condition of the multigrid preconditioner in order to satisfy requirements of a 
preconditioner of the PCG method. Next numerical experiments show a behavior of the MGCG 
method and that the MGCG method is superior to both the ICCG method and the multigrid 
method in point of fast convergence and high parallelism. This fast convergence is understood in 
terms of the eigenvalue analysis of the preconditioned matrix. From this observation of the 
multigrid preconditioner, it is realized that the MGCG method converges in very few iterations and 
the multigrid preconditioner is a desirable preconditioner of the conjugate gradient method. 

1 INTRODUCTION 

The typical numerical methods of a king-size system of linear equations, after discretization of 
the partial differential equations, are the preconditioned conjugate gradient method (PCG method) 
and the multigrid method [12]. The conjugate gradient method is valued in that it suits to parallel 
computing and even ill-conditioned problems can be easily solved with the help of a good 
preconditioning. 

This paper considers an efficient preconditioner and proposes a multigrid preconditioned 
conjugate gradient method (MGCG method) which is the conjugate gradient method with the 
multigrid method as a preconditioner. The combination of the multigrid method and the conjugate 
gradient method was already considered. Kettler and Meijerink [7] and Kettler [8] treated the 
multigrid method as a preconditioner of the conjugate gradient method. However this paper 
formulates the MGCG method more generally than these and requirements of the multigrid 
preconditioner are studied. On the other hand, Bank and Douglas [2] treated the conjugate gradient 
method as a relaxation method of the multigrid method. Braess [3] considered these two 
combinations and reported that the conjugate gradient method with a multigrid preconditioning is 
effective for elasticity problems. 
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We study requirements of the valid multigrid preconditioner and evaluate this preconditioner by 
some numerical experiments and eigenvalue analysis. Especially, eigenvalue analysis is a more direct 
and more reasonable criterion than convergence rate, since the number of iterations of the conjugate 
gradient method until convergence depends on the eigenvalues’ distribution of the preconditioned 
matrix. In Sections 2 and 3, the preconditioned conjugate gradient method and the multigrid 
method which axe the basis of this paper are briefly explained. Section 4 discusses the requirements 
of the valid two-grid preconditioner for the conjugate gradient method; then in Section 5, it is 
extended to the requirements of the multigrid preconditioner. In Section 7, numerical experiments 
show that the MGCG method converges with very few iterations even for ill-conditioned problems. 
In Section [8], eigenvalue analysis is performed, and it is realized why the MGCG method can easily 
solve the problem that the ordinary multigrid method itself does not converge rapidly. When the 
multigrid method is used as a preconditioner of the conjugate gradient method, it becomes quite an 
effective and desirable preconditioner of the conjugate gradient method. 

2 THE PRECONDITIONED CONJUGATE GRADIENT METHOD 


If a real n x n matrix A is symmetric and positive definite, the solution of a linear system 
Ax = f is equivalent to minimization of the quadratic function 

Q{x) = ^ x t Ax - f T x. ( 1 ) 

The conjugate gradient method is one of the minimization methods and uses A-conjugate vectors as 
direction vectors which are generated sequentially. Theoretically this method has the striking 
property that the number of steps until convergence is at most n steps. This method can be 
adapted successfully to the parallel and vector computation, since one CG iteration requires only 
one product of the matrix and the vector, two inner products, tree linked triads, two scalar divides 
and one scalar compare operation. 

Next the preconditioned conjugate gradient method is explained. Let U be a nonsingular matrix 
and define A = U AU T ; then solve Ax = f using plain conjugate gradient method. Let x° be an 
initial approximate vector; then an initial residual r° is r*° = / — Ax 0 . Let M = U T U, f° = Mr 0 
and an initial direction vector p° = f°. The PCG algorithm is described by Program 1. 

The matrix M is a precondition matrix and this paper focuses on this computation. A new 
proposal is the PCG method exploiting the multigrid method as a preconditioner. 

On the other hand, the matrix M should satisfy some conditions: symmetric and positive 
definite. Therefore if the matrix of the multigrid method is symmetric and positive definite, it is 
reasonable to use the multigrid method as a preconditioner of the CG method. In Sections 4 and 5, 
the conditions of the multigrid preconditioner in order to satisfy the requirements of a 
preconditioner of the conjugate gradient method are investigated. 
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i = 0; 

while ( ! convergence ) { 

«» = {r^T^KPiyApi); 
x i+ i = Xi + aiPi] 

’’i+i = I'i ~ (XiAPii 
convergence test; 

r, + i = Mri+i ; // preconditioning 

/?» = (fj + i,r j+ i)/(fj,ri); 

Pt+I = ^t+1 d" fiiPii 

i++; 

} 


Program 1. The PCG iteration 
3 THE MULTIGRID METHOD 


In the iterative methods, the frequency components of the residual are reduced most rapidly on 
the grid corresponding to them. The multigrid method makes good use of this characteristic and 
exploits a lot of grids to converge as rapid as possible. 

These grids are leveled and numbered from the coarsest grid. This number is called the level 
number. If the multigrid method is applied to the solver of linear equations, the residual is reduced 
moving it from grid to grid. The basic element of the multigrid method is the defect correction 
principle. The defect correction scheme consists of three processes: pre-smoothing process, coarse 
grid correction and post-smoothing process. In the smoothing process, various methods, such as 
ILU, ADI and zebra relaxation, are proposed. One purpose of this research is, however, formation of 
an efficient method with high parallelism. Thus an iterative method with high parallelism, such as 
the damped Jacobi method or a multi-color symmetric SOR method (SSOR method), is used as the 
smoothing method. 

An operation of transferring a vector on a finer grid to a vector on a coarser grid is called 
restriction , and an opposite operator is called prolongation. A matrix presenting the operation of 
restriction is written r in this paper, and prolongation is p. 

In the following section, the equation of grid level i is described as 

LiXi = fi 

and restriction is defined by adjoint of prolongation. That is, 

r = bp T , 


where 6, a scalar constant, is satisfied. 


623 




4 THE TWO-GRID PRECONDITIONER 


This section and the next section examine whether the multigrid method suits a preconditioner 
of the PCG method. First it is shown that two kinds of two-grid methods, one with pre-smoothing 
and no post-smoothing and the other with both pre-smoothing and post-smoothing, satisfy the 
conditions of a preconditioner: the matrix of the two-grid method is symmetric and positive 
definite. Next it is shown that V-cycIe and W-cycle multigrid methods also hold. 

A linear equation, = f ,, is concerned. If R is a matrix of a relaxant calculation and u is an 
approximate vector, one two-grid iteration can be shown by matrix form in Table 1. 


u = H m u + Rf 
d = r ( L t u - f ) 
v = Lj\d 
u — u — pv 
u — H m u + Rf 


// pre- smoothing 
// coarse grid correction 


// post-smoothing 


Table 1. The two-grid iteration 


In this paper the relaxant calculation is an iterative method with high parallelism, and the 
matrix R is defined as follows. Let Li be an n x n nonsingular, symmetric matrix and be split as 

L, = P-Q, ( 2 ) 

where P is a nonsingular matrix and the symmetric part of P -|- Q is positive definite. For example, 
in the case of the point Jacobi method, Pis a diagonal matrix containing diagonal elements of L\. 
Then the i’th approximate vector u' is updated such as 

u i+1 = P~ l Qu i + P~ l f. (3) 

If an initial approximate vector U® is zero vector and w. iterations are done, R is equal to 

m- i 

R = 53 R~ l i R = P~ l Q • (4) 

i= 0 

H is called an iterative matrix. 


4.1 The two-grid preconditioner with pre-smoothing only 

First consider a no post-smoothing case. The matrix of one iteration of Table 1 equals 

M = (I-pLT_\rL,)R+pL;\r 

= R + pLi} ir (I-L,R). (5) 

Then the following theorem holds. 
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Theorem 1 The matrix Lff x is symmetric and positive definite, and N = I — LiR. If the matrix N 
and P are symmetric, the matrix M of Eq. (5) is symmetric in the N-energy inner product. If the 
matrix N is symmetric and nonsingular, the matrix P is symmetric and m is even; then the matrix 
M is positive definite in the N-energy inner product, provided that N-energy inner product 
{x,v)n = (x,Ny). 


Proof. Since N is symmetric, (I — L X R) T = I — LiR. Therefore 

I - R t L, = I - L,R. 

Since P is symmetric, the matrix R is also symmetric. Then 

RL, = UR. (6) 

And 

(x,My)jy = x T NRy + x T NpLj} x r{I — L x R)y 

= x T {I-UR)Ry + x T {I-UR)pLj} x r{I-UR)y. (7) 


Besides 


(Mx,y)x — x T M T Ny 

= x T RNy + x T (I — LiR)pLj} ] rNy 

= x T (I - L t R)Ry + x T (I - LiR)pLj\r(I - L t R)y (because of Eq. (6)) 

= ( x,My) N . (8) 

Therefore the matrix M is symmetric in the iV-energy inner product. 


Next, it is shown that the matrix M is the positive definite in the iV-energy inner product. It is 
equal to (x,Mx)n > 0. Then 

N = I -UR 
= I-RU 

m—i 

= /- Y.( p ~'Q) ip -'( p -Q) 

t=0 

= (. P~ l Q) m 
= H m . 


Thus 


NM = (I - L,R){R + pLjl x r (I - UR)} 

= H m R + H m pLjf x rH m . 

Since P is symmetric and nonsingular and L\ is symmetric and positive definite, then 
H = P~ l Q = I — P~ l Li has real eigenvalues. Hence if m is even, H m is positive definite. If P + Q 
is positive definite and m is even, then R is positive definite (see [11]). Therefore H m R is positive 
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definite. Since H m is symmetric and pLjfpr is semi-positive definite, H m pL l f l rH m is semi-positive 
definite. Thus NM is positive definite. 0 

The iterative method which holds the assumption of Theorem 1 is the damped Jacobi method. 
From this theorem, the two-grid preconditioner with the damped Jacobi method as a relaxant 
calculation fills the conditions of the preconditioner of the CG method which uses the iV-energy 
inner product instead of the usual inner product. 

4.2 The two-grid preconditioner with both pre-smoothing and post-smoothing 

Next consider the two-grid iteration with both pre-smoothing and post-smoothing. Suppose the 
pre-smoothing and the post-smoothing are the same method. Then the matrix of one two-grid 
iteration in Table 1 equals 

M = H m {(I — pLj} x rL{)R + pLj\r} + R 

= H m R + R + H m pLjl x r {I - L,R). (9) 

However since P and Q axe symmetric, 

I-L,R={QP~ l ) m = {H T ) m . 

Therefore the matrix M of Eq. (9) is rewritten as 

M = H m R + R + H m pL(} l r(H T ) m . (10) 

Then the following theorem is satisfied. 

Theorem 2 The matrix L(\ is symmetric and positive definite. If the matrix P is symmetric, the 
matrix M of Eq. (10) is symmetric and positive definite. 

Proof. Since the matrix P is symmetric, the matrix R is also symmetric. Thus 

M t = R(H T ) m + R + H m pLj\r ( H T ) m . 

Now 

m — 1 

H m R = H m Y, WP- 1 . 

»=0 

771-1 

R(H T ) m = Yj H i P~ 1 (H T ) m - 

i=0 

Moreover since P is symmetric and H = P~ l Q , then P~ l H T = HP~ l . Therefore 

H m R= R(H T ) m . 
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After all, the matrix M is symmetric. Next show that the matrix M is positive definite. 

M = H m R + R+H m pLjl x r{H T ) m 

2m— l 

= £ WP- 1 + H m pLf} x r ( H T ) m . (11) 

t=0 

2m — l 

Since the first term of right hand expression ^ H'P~ l of Eq. (11) is the matrix after 2m times 

i=0 

iteration, it is positive definite if P + Q is positive definite. Since Lj\ is positive definite, 

H m pLj} x r (H T r is semi-positive definite. Therefore M is positive definite. □ 

The iterative methods which hold the assumption of Theorem 2 are the damped Jacobi method, 
Red-Black Symmetric Gauss-Seidel method (RB-SGS method), multi-color SSOR method, ADI 
method and so on. From this theorem, the two-grid preconditioner with one of these iterative 
methods as a relaxant calculation fulfills the conditions of the preconditioner of the CG method. 

5 THE MULTIGRID PRECONDITIONER 


In the previous section the possibility of two kinds of two-grid preconditioners is considered. In 
the following, only the latter two-grid preconditioner, with both pre-smoothing and post-smoothing, 
is discussed. However the same discussion can be applied to the former two- grid preconditioner. In 
this section, extension to the multigrid preconditioner is argued. The following theorem holds. 

Theorem 3 If assumptions of Theorem 1 and 2 are satisfied, all MG(m , n) methods (m, n > 1) 
satisfy conditions of a preconditioner of the CG method, where m is a multigrid cycle and n is the 
number of iterations of the smoothing method. 


Proof. The matrix Mi of the V-cycle multigrid method can be defined as 
Mo = Lq 1 or R* 

Mi = H m Ri + R i + H m pM i . x r{H T ) m . {i > 1) 

M 0 is symmetric and positive definite. If M, is symmetric and positive definite, M i+X is also 
symmetric and positive definite because of Theorem 2. By mathematical induction, every 
Mi{i > 0) is symmetric and positive definite. Therefore the V-cycle multigrid method can be used 
as a preconditioner. 

Next the W-cycle multigrid method is considered. If the matrix iVj”^ is the multigrid method 
with n recursive calls of the multigrid method on level number /— 1 as the solution on the coarse 
grid, Nj^ is defined as 

N{ n) = L^ovRo 

N ( i n) = H' mi {H m Ri + Ri + H m pN^\r (H T ) m }, (i > 1) 

i=0 
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where = H 2m — H m pN-"\ rL{H m . is symmetric and positive definite. If IV-"} is symmetric 

and positive definite, H m Ri + Ri + H m pN^\r (H T ) m is symmetric and positive definite by 
Theorem 2. Thus N is also symmetric. And because of p(H mg ) < 1 by [?]> N < T is positive 
definite. The W-cycle multigrid method is the case of n = 2. Therefore the W-cycle multigrid 
method and all MG(n, m) (m,n > 1) satisfy the conditions of the preconditioner. □ 

6 THE MGCG METHOD 


In the previous section, the multigrid preconditioner which is valid for a preconditioner of the 
CG method is considered. When only pre-smoothing is performed, the multigrid preconditioner 
with an even number of iterations of the damped Jacobi smoothing can become a preconditioner of 
the conjugate gradient method with the iV-energy inner product instead of the usual inner product. 
When both pre-smoothing and post-smoothing are performed, the multigrid preconditioner with 
RB-SSOR smoothing, ADI method and so on, fulfills the requirements of a preconditioner of the 
conjugate gradient method. Thus the multigrid preconditioned conjugate gradient method (MGCG 
method) is mathematically valid. There are several variations of this preconditioner. If m is a cycle 
of the multigrid method, / is a relaxant method, n is the number of iterations of the relaxant 
method and g is the number of grids, the MGCG method is specified as MGCG(/, m,n,g). But g is 
an optional parameter and if this parameter is omitted, all available grids are used. For example, 
MGCG(jRJ 9, 1, 2) is the MGCG method of the V-cycle multigrid preconditioner with two iterations 
of the Red-Black SSOR smoothing. 

7 NUMERICAL EXPERIMENTS 
7.1 Problems 


A two-dimensional Poisson equation with Dirichlet boundary condition: 

—V (fcVu) = / in D = [0, 1] x [0, 1] 
with u — g on dft, 

where A: is a real function, is considered. The equation is defined by a diffusion constant k, a source 
term / and a boundary condition g. Numerical experiments are performed in the following two 
conditions. 


Problem 1 Diffusion constant is uniform and source term is equal to 0. Boundary condition is 
<7 = 0 except y = 1 and g = 3z(l — x) on y = 1. 

Problem 2 Diffusion constant and source term are depicted by Figs. 1 and 2. Boundary condition 
g is always equal to 0. 
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Figure 1. Diffusion constant of problem 2 Figure 2. Source term of problem 2 



Problem 1 is a simple case, and the multigrid method is expected to converge efficiently. The 
multigrid preconditioner is also expected to be efficient. Problem 2 has a non-uniform diffusion 
constant and the area with a large diffusion constant looks like a letter ‘T’; therefore it has a rich 
distribution of eigenvalues of the problem matrix, which is investigated in the next section. 
Moreover since a source term is complex, it does not happen that specific iterative methods, such as 
ICCG method and MICCG method, accidentally converge very rapidly. 

These problems are discretized to three kinds of meshes: 64 x 64, 128 x 128 and 256 x 256, by 
the finite element method. These coefficient matrices become symmetric, positive definite and block 
tridiagonal. 


7.2 Solutions 


In numerical experiments, three methods: the MGCG(il.B, 1, 2) method, the ICCG(1, 2) 
method and the MG(1, 2) method, are compared. The ICCG(1, 2) method is the PCG method 
with the incomplete Cholesky decomposition having an additional one line to the original problem 
sparse matrix. The MG(1, 2) method is the identical method to the multigrid preconditioner of the 
MGCG(RB, 1, 2) method. 

Numerical experiments are performed on the HP9000/720 and the program is written by C++ 
with original vector and matrix classes. 
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7.3 Convergence of the MGCG method 



MGCG(i?B, 1,2) 


ICCG(1, 2) 

MG 

(1,2) 

size 

# of iter. 




iter. 

time 

iter. 

time 


5 


4 


38 

1.19 

H 

0.65 

■WEI 

5 

3.16 

5 

4.58 

72 

10.88 


4.05 

255 2 

5 

15.8 

5 

23.7 

134 

89.5 

H 

20.2 


(HP9000/720; C++) 


Table 2. Problem 1 



MGCG(i?B, 1, 2) 


ICCG (1,2) 

MG(1, 2) 

size 

# of iter. 

time(sec.) 


time 

iter. 

time 


time 

63 2 

9 


8 

■ 

53 

1.65 


13.4 

127 2 

9 

5.54 

8 


103 

15.49 

B 

75.3 

255 2 

9 

27.8 

8 

37.4 

200 

133.0 

122 

341.5 


(HP9000/720; C++) 


Table 3. Problem 2 

Tables 2 and 3 are results of these numerical experiments. The number of it erat ions and the time 
of each method until convergence are measured. The number of iterations of the MGCG method 
and the ICCG method is that of CG iterations and the number of iterations of the multigrid method 
is that of V-cycle iterations. From results of the two problems, the following points are notable: 

• The MGCG method converges with very few iterations. 

• The number of iterations of the MGCG method does not increase when a mesh size is larger. 

• Even for complex problems, such as problem 2, the MGCG method converges fast. 

The first item is discussed by an eigenvalue analysis in the next section. From the second item, the 
MGCG method is advantageous over the ICCG method even as large as the mesh size is. It is a 
principle of the multigrid method that the number of iterations does not depend upon the mesh 
size. If the problem is simple such as problem 1, the multigrid method converges very fast; however, 
in complex problems, such as problem 2, it converges very slowly. To avoid this, the multigricfcXI .• ' 
method should have the stronger relaxation method, but the stronger relaxation method has poor 
parallelism. Moreover in problem 2, it is considered that the locking effect [?] has occurred. From 
the third item, the MGCG method is also superior to the multigrid method as a result of stably fast 
convergence and high parallelism. 
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8 EIGENVALUE ANALYSIS 


In order to study the efficiency of the multigrid preconditioner, the eigenvalue distribution of a 
coefficient matrix after preconditioning is examined. The number of iterations of the conjugate 
gradient method until convergence depends upon an initial vector, a distribution of eigenvalues of a 
coefficient matrix and a right-hand term, but due to a good initial vector and a simple right-hand 
term, the conjugate gradient method happens to converge fast unreasonably, so the eigenvalue 
distribution is investigated. The problem is the same problem in Section 7 and the area is 
discretized to the mesh of 16 x 16 by the finite element method. The condition number of this 
coefficient matrix is 5282.6. 

A matrix after the multigrid preconditioning is calculated as follows. The matrix M of Eq. (5) 
or (10) is Cholesky decomposed as M = U T U, then eigenvalues of the matrix U LiU t is investigated. 
On the other hand the matrix using the ICCG method is calculated as follows. The matrix L\ is 
incomplete Cholesky decomposed as L t = S T S - T, and the general eigenvalue problem 
L(X = XS T Sx is solved in order to examine eigenvalues after preconditioning. 



Figure 3. Eigenvalue distribution of a problem 
matrix 



Figure 4. Eigenvalue distribution after 
preconditioning 


The eigenvalue distribution of the problem matrix is shown in Fig. ??. The horizontal x axis is 
the order of the eigenvalues and the vertical y axis values are the eigenvalues. This vertical axis is in 
a log scale. The eigenvalue distribution of the matrix after preconditioning is shown in Fig. ??. This 
vertical axis is in a linear scale. In order to compare, preconditioning is carried out in both the 
multigrid method and the incomplete Cholesky decomposition. 

The eigenvalue distribution of the multigrid preconditioner is effective for the conjugate gradient 
method as the following points: 
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1. Almost all eigenvalues are clustered around 1 and a few eigenvalues are scattered between 1 
and 0. 

2. The smallest eigenvalue is larger than the ICCG method. 

3. Condition number is decreased. 

The first item is no problem for the conjugate gradient method. All of these characteristics are 
desirable to accelerate the convergence of the conjugate gradient method. In problem 1, there are 
no scattered eigenvalues. So the multigrid method converges efficiently, however in problem 2, the 
scattered eigenvalues prevent the ordinary multigrid method from converging rapidly. Therefore 
using the multigrid method as a preconditioner of the conjugate gradient method is quite important. 

9 CONCLUSION 

This paper investigates the conjugate gradient method with a multigrid preconditioner (MGCG 
method). Necessary conditions of a preconditioning matrix of the conjugate gradient method are 
symmetric and positive definite. First two kinds of two-grid preconditioners are considered and 
conditions of both preconditioners are given in order to satisfy necessary conditions of a 
preconditioner. Secondly extension to the multigrid preconditioner is carried out and conditions for 
a valid multigrid preconditioner are also given. Thirdly numerical experiments are performed and 
the MGCG method has faster convergence and a more effective method than both the ICCG 
method and the multigrid method. Finally eigenvalue analysis is performed in order to verify the 
effect of the multigrid preconditioner. It concludes that the multigrid preconditioner is an excellent 
preconditioner and it improves the number of the CG iterations remarkably. Consequently the 
MGCG method has the following properties: 

• The number of iterations does not increase even when a mesh is finer. 

• Even in the case that the problem is ill-conditioned, the MGCG method is effective. 

• The distribution of the eigenvalues of the matrix after preconditioning is suited to the 
conjugate gradient method. 

• The MGCG method has high parallelism. 


The multigrid method roughly solves any problems, since almost all eigenvalues of Section ?? are 
clustered around the unity, but a few scattered eigenvalues prevent fast convergence. The conjugate 
gradient method hides the defect of the multigrid method. Therefore the MGCG method becomes 
an efficient method. Parallelization of the MGCG method and implementation on the 
multicomputers are beyond the scope of this paper, so this facility is no more mentioned. However 
since the MGCG method has high parallelism and fast convergence, this method is a very promising 
method as the solution of a large-scaled sparse, symmetric and positive definite matrix. 
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